05 Machine Learning ¶
Right click to download this notebook from GitHub.
With the data preparation complete, this step will demonstrate how you can configure a scikit-learn or dask_ml pipeline, but any library, algorithm, or simulator could be used at this stage if it can accept array data. In the next step of the tutorial, Data Visualization you will learn how to visualize the output of this pipeline and diagnose as well as ensure that the inputs to the pipeline have the expected structure.
Recap: Loading data ¶
import intake
cat = intake.open_catalog('../catalog.yml')
l5_da = cat.l8.read_chunked()
Subsetting data ¶
To speed things up for tutorial purposes, we'll use a subset of these data in the following examples. There are many ways to subset data in xarray . Here we select the central third of the data using index selection, which will be around 1500 x 1500 pixels.
nbands, ny, nx = l5_da.shape
bounds = int(2*ny/5), int(3*ny/5), int(2*nx/5), int(3*nx/5)
bounds
l5_da = l5_da[:, bounds[0]:bounds[1], bounds[2]:bounds[3]]
l5_da
import hvplot.xarray
l5_da.hvplot(kind='image', x='x', y='y', groupby='band', datashade=True, width=450, height=400, cmap='greys')
Reshaping Data ¶
We'll need to reshape the image to be how dask-ml / scikit-learn expect it:
(n_samples, n_features)
where n_features is the number of bands and n_samples is the total number of pixels in each band. Essentially, we'll be creating a bag of pixels out of each image, where each pixel has multiple features (bands), but the ordering of the pixels is no longer relevant. In this case we start with an array that is n_bands by n_y by n_x (7, 2000, 2000) and we need to reshape to an array that is
(n_samples, n_features)
(4e6, 7). We'll first look at using NumPy, then Xarray.
Numpy ¶
Data can be reshaped at the lowest level using NumPy, by getting the underlying values from the
xarray.DataArray
, and using flatten and transpose to get the right shape.
import numpy as np
arr = l5_da.values
arr.shape
Since we want to flatten along the x and y but not along the band axis, we need to iterate over each band and flatten the data.
flattened = np.array([arr[i].flatten() for i in range(arr.shape[0])])
flattened
flattened.shape
We can reorder the dimensions using
.transpose
sample_by_feature = flattened.transpose()
sample_by_feature
sample_by_feature.shape
Since
numpy.array
s are not labeled data, the semantics of the data are lost over the course of these operations, as the necessary metadata does not exist at the NumPy level.
Xarray ¶
By using
xarray
methods to flatten the data, we can keep track of the coordinate labels ('x' and 'y') along the way. This means that we have the ability to reshape back to our original array at any time with no information loss.
flattened_by_band = l5_da.stack(z=('x','y'))
flattened_by_band
We can reorder the dimensions using
Dataset.transpose
:
sample_by_feature = flattened_by_band.transpose('z', 'band')
sample_by_feature
Now we have the data in the shape that we are looking for: a long array of pixels for each band. As a sanity check we can take a look at the plain
np.array
:
X = sample_by_feature.values
X
np.expand_dims(X, 2).shape
sample_by_feature.expand_dims(dim='e', axis=2)
# Exercise: Try removing the extra axis using np.squeeze or .squeeze on the xarray object
Rescale: ¶
Rescale (standardize) the data to input to the algorithm since the ML pipeline that we have selected expects input values to be small.
(X - X.mean()) / X.std()
rescaled = (sample_by_feature - sample_by_feature.mean()) / sample_by_feature.std()
rescaled.compute()
NOTE:
Since the the xarray object is in dask, the actual computation isn't performed until
.compute()
is called.
# Exercise: Inspect the numpy array at rescaled.values to check that it matches the numpy array above.
ML pipeline ¶
The Machine Learning pipeline shown below is just for the purpose of understanding the shaping/reshaping of the data. In practice you will likely be using a more sophisticated pipeline. Here we will use a version of SpectralClustering from dask_ml that is a scalable equivalent to operations from Scikit-learn that cluster pixels based on similarity (across all bands, which makes it spectral clustering by spectra!).
from dask_ml.cluster import SpectralClustering
from dask.distributed import Client
client = Client(processes=False)
client
Now we will compute and persist the rescaled data to feed into the ML pipeline. Notice that X has the shape:
n_samples, n_features
as discussed above.
X = client.persist(rescaled)
X.shape
clf = SpectralClustering(n_clusters=4, random_state=0, gamma=None,
kmeans_params={'init_max_iter': 5},
persist_embedding=True)
%time clf.fit(X)
# Exercise: Open the dask status dashboard and watch the workers in progress.
labels = clf.assign_labels_.labels_.compute()
labels.shape
Un-flattening ¶
Once the computation is done, the output can be used to create a new array with the same structure as the input array. This new output array will have the coordinates needed to be unstacked similarly to how they were stacked. One of the main benefits of using
xarray
for this stacking and unstacking is that allows
xarray
to keep track of the coordinate information for us.
template = sample_by_feature[:, 0]
template
NOTE: Since the original array is n_samples by n_features (4_000_000, 6) and the output only contains one feature (4_000_000,), the template structure for this data needs to have the shape (n_samples). We achieve this by just taking one of the bands.
output_array = template.copy(data=labels)
output_array
With this new output array in hand, we can unstack back to the original dimensions:
unstacked = output_array.unstack()
unstacked
l5_da.sel(band=4).hvplot(x='x', y='y', datashade=True, cmap='greys', width=450, height=400).relabel('Image') + \
unstacked.hvplot(x='x', y='y', datashade=True, cmap='Category10', width=450, height=400).relabel('Clustered')
Geographic plot ¶
The plot above is useful and quick to generate, but it isn't referenced against the underlying geographic coordinates, which is crucial if we want to overlay the data on any other geographic data sources. Adding the coordinate reference system in the hvplot method, ensures that the data is properly positioned in space. This geo-referencing is made very straightforward because of the way
xarray
persists metadata. We can even add tiles underneath.
import geoviews.tile_sources as gts
gts.ESRI * unstacked.hvplot(x='x', y='y', datashade=True, geo=True, height=500, cmap='Category10')
# Exercise: Try adding a different set of map tiles. Use tab completion to find others.
Next: ¶
Now that your analysis is complete, you are ready for some more information about Data Visualization you will learn how to visualize the output of this pipeline and diagnose as well as ensure that the inputs to the pipeline have the expected structure.